L’énergie éolienne constitue un moyen propre et renouvelable pour produire de l’électricité, utilisée depuis plusieurs siècles (moulins à vent), cette dernière s’est particulièrement développée depuis l’avènement de l’éolien électrique (seconde moitié du XXe siècle).
Premier acteur éolien en France, ENGIE fait du développement de l’éolien une de ses priorités et vise un objectif de 2 000 MW de puissance installée.
Détecter des écarts anormaux entre la production électrique attendue et la production réalisée aide à mieux optimiser la production de chaque éolienne. Dans l’optique d’améliorer la production d’électricité éolienne nous allons analyser les données fournies par ENGIE sur la plateforme Data challenge dans le cadre du cours ”Analyse de données et apprentissage (M2 Ingénierie Statistique et Data Science UPMC&ISUP)”.
Dans un premier temps nous étudierons les données, puis nous chercherons à mettre en place un premier modèle, nous comparerons ensuite ce dernier avec d’autres modèles et nous essayerons d’améliorer le modèle sélectionné.
#Descriptif des données
Les données qui seront utilisées pour ce projet sont composées de 617386 observations réparties en deux fichiers. Le premier input contient les 78 variables explicatives et le deuxième fichier contient la variable à expliquer Target. Chaque ligne de notre jeu de données est définie par un identifiant unique ID, liée à une éolienne MAC_CODE et un pas de temps Data_time.
L’objectif, est de prédire à partir de l’ensemble de données la puissance active des éoliennes.
#Pre-traitement des données
Nous commençons par une brève description des données via les fonctions summary, glimpse, etc. Nous remarquons que la variable MAC_CODE est une variable factorielle, le reste des variables sont quantitatives continues. La variable Target étant continue, nous cherchons donc à faire une régression.
Le jeu de données contient 469944 valeurs manquantes réparties sur 109416 lignes. Il existe plusieurs méthodes pour l’imputation des données manquantes: plus proche voisin médian, Amélia, MissForest …etc.
Comme le jeu de données est volumineux et la capacité de notre machine est limitée, nous allons tout simplement supprimer toutes les lignes contenant des NA’s.
Nous avons remarqué que pour toutes les variables nous avons le min (minimum),le max (maximum) et le std (ecart-type). D’abord, nous allons prendre quelques exemples de variables: Pitch_angle, Hub_temperature, Nacelle_angle et Grid_frequency. Puis, visualiser leurs lien avec leurs min, max et std respectifs.
D’après la matrice de corrélation ci-dessus, les variables Pitch_angle, Hub_temperature, Nacelle_angle et Grid_frequency sont très fortement corrélées avec leurs min et max respectifs. En revanche, ce n’est pas le cas avec leurs std.
Maintenant, prenons l’exemple de la variable Hub_temperature et compraons sa distribution à celles de Hub_temperature_min, Hub_temperature_max et Hub_temperature_std.
D’après les graphiques ci_dessus, les variables Hub_temperature,Hub_temperature_min et Hub_temperature_max ont une même distribution de densité et prennent leurs valeurs dans l’intervalle \([0,50]\). Par contre Hub_temperature_std a une distribution différente des précédente et elle prend ses valeurs dans l’intervalle \([0,0.5]\). Afin d’éviter la redondance, il suffit de prendre Hub_temperature et Hub_temperature_std.
Nous avons eu des résultats sémilaire pour toutes les autres variable, donc pour la suite, nous allons supprimer toutes les variables min et max.
Maintenant, observons la matrice de corrélation des variables restantes après suppression des variables min et max.
D’après la matrice de corrélation:
\(\bullet\) Les variables Generator_converter_speed etGenerator_speed sont très fortement corrélées.
\(\bullet\) Les variables Generator_converter_speed_std etGenerator_speed_std sont très fortement corrélées.
\(\bullet\) Les variables Nacelle_angle_c etNacelle_angle sont très fortement corrélées.
Nous allons donc supprimer les variables: Generator_converter_speed,Generator_converter_speed_std et Nacelle_angle_c.
L’objectif de cette partie est d’étudier tous les points atypiques et éventuellement supprimer des individus.
D’après le Normal Q-Q plot, un nombre de points important ne sont pas alignés sur la première bissectrice de plus ils dépassent les bornes \([-2,2]\). On ne peut pas dire que tous ces points sont aberrants et se permettre de les supprimer.
MAC_CODED’aprés l’histogramme de la variable MAC_CODE qui désigne l’éolienne sur laquelle les données ont été mesurés, on constate que les proprtions des quatres éolienne sont prêsque les mêmes avec une légère hausse pour l’éolienne WT3.
TargetD’après l’histogramme ci-dessus, la variable TARGET prend ses valeurs entre \(-20\) et \(2000\), avec une concentration autour de zéro.
Les graphiques ci-dessus représentent les histogrammes des variables Rotor_speed, Hub_temperature et Generator_speed_std.
D’abord la variable Rotor_speed ne prend que des valeurs positives, dont plus de 80% sont dans [0,2.5[ \(\bigcup\) [8,16]. Ensuite, nous remarquons que Hub_temperature est unimodale et asymétrique.
Enfin, la variable Generator_speed_std prend que des valeurs positives, de plus elle est unimodale et asymétrique.
D’après les boxplot et les graphes de densités ci-dessus, nous remarquons que les quantiles de la variable TARGET sont légèrement différents d’une éolienne à l’autre, de même pour les densités. La variable TARGET n’a donc pas la même loi pour les quatres éoliennes.
| Variable | Disperssion | ||
|---|---|---|---|
| Rotor_speed_std | 7.550396e-01 | ||
| Grid_voltage_std | 4.583308e+00 | ||
| Generator_speed_std | 7.925428e+01 | ||
| Date_time | 4.881169e+04 |
Tableau illustrant la disperssion de quelques variables
Après avoir observé la disperssion des variables, nous avons constaté que les échelles sont très différentes. Nous avons des variables de l’ordre \(10^{-1}\) jusqu’à l’ordre \(10^{4}\), il faut donc normaliser.
Afin d’identifier les variables qui expliquent le mieux notre jeux de données,les variables corrélées et réduire les dimensions, nous allons appliquer une ACP normée. Pour le calcul de l’ACP nous allons utiliser les packages “FactoMiner” et “factoextra”.
Les valeurs propre mesurent la quantité de variance expliquée par chaque axe principal, nous allons donc les examiner afin de déterminer le nombre de composantes principales à prendre en compte.
| Composante | Valeur propre | Pourcentage de variance | Pourcentage cumulé de variance | ||
|---|---|---|---|---|---|
| comp 1 | 9.785503e+00 | 2.878089e+01 | 28.78089 | ||
| comp 2 | 3.958601e+00 | 1.164294e+01 | 40.42383 | ||
| comp 3 | 3.035529e+00 | 8.928027e+00 | 49.35186 | ||
| comp 4 | 2.241573e+00 | 6.592862e+00 | 55.94472 | ||
| comp 5 | 1.518758e+00 | 4.466934e+00 | 60.41166 | ||
| comp 6 | 1.501787e+00 | 4.417020e+00 | 64.82868 | ||
| comp 7 | 1.105757e+00 | 3.252226e+00 | 68.08090 | ||
| comp 8 | 1.048219e+00 | 3.082998e+00 | 71.16390 | ||
| comp 9 | 1.000882e+00 | 2.943771e+00 | 74.10767 | ||
| comp 10 | 9.868843e-01 | 2.902601e+00 | 77.01027 |
Tableau des valeurs propres de quelques composantes
Dans l’ACP normée, la règle pour choisir la dimension est:
\[q=max\left\{l/\lambda_l \geq 1\right\}\] D’après le tableau ci-dessus, le nombre d’axes principeaux à conserver est: \(q=9\).
Une autre manière de déterminer le nombre de composantes principales est de regarder le graphique des valeurs propres (scree plot).
D’après le graphique à partir de \(q=9\) ou \(q=10\) les valeurs propres restantes sont toutes relativement petites, on peut donc choisir de garder les dix premiers axes.
\(\bullet\) Nuage des variables
Nous utiliserons le cercle de corrélation pour visualiser les variables:
Le graphique ci-dessus illustre les relations entre toutes les variables, prenons par exemple:
\(\circ\) Les variables Gearbox_bearing_2_temperature_std,Generator_bearing_2_temperature_std et Grid_voltage sont regroupées, elles sont donc postivement corrélées. En revanche, elles sont négativement corrélées avec la variable Pitch_angle.
\(\circ\) Les variables Hub_temperature et Gearbox_oil_sump_temperature_std ne sont pas corrélées.
\(\circ\) Les variables Hub_temperature et Generator_bearing_2_temperature sont loin de l’origine, elles sont donc bien représentées par l’ACP.
\(\bullet\) Qualité de représentation des variables
La qualité de représentation des variables est mesurée par \(cos^2\). Ci-dessus le bar-plot correspondant:
Les variables Generator_bearing_1_temperature,Gearbox_oil_sump_temperature et Nacella_temperature ont un \(cos^2\) elevé sont bien représentées par les axes principaux 1-2. En revanche, les variables Grid_voltage, Nacella_angle et Grid_frequency ont un \(cos^2\) très faible, elles sont donc mal représentées par les axes principaux 1-2.
\(\bullet\) Contribution des variables
| Dim.1$quanti | Corrélation | |
|---|---|---|
| Generator_speed_std | 0.76032684 | |
| Rotor_speed_std | 0.75877857 | |
| Gearbox_bearing_1_temperature_std | 0.62360100 | |
| Generator_stator_temperature | -0.087450697 | |
| Gearbox_oil_sump_temperature | -0.106711545 | |
| Nacelle_angle_std | -0.113520239 |
| Dim.2$quanti | Corrélation | |
|---|---|---|
| Gearbox_oil_sump_temperature_std | 0.524247615 | |
| Rotor_speed | 0.468734326 | |
| Generator_speed | 0.468601377 | |
| Pitch_angle_std | -0.136354131 | |
| Pitch_angle | -0.203310594 | |
| Gearbox_inlet_temperature | -0.206463311 |
| Dim.3$quanti | Corrélation | |
|---|---|---|
| Generator_speed_std | 0.76032684 | |
| Rotor_speed_std | 0.75877857 | |
| Gearbox_bearing_1_temperature_std | 0.62360100 | |
| Gearbox_inlet_temperature | -0.10513183 | |
| Grid_voltage | -0.10792038 | |
| Grid_frequency | -0.10919804 |
Tableau illustrant une partie de la sortie `dimdesc`
En visualisant les sorties de dimdesc, nous observons:
Plusieurs variables sont très liées positivement avec l’axe 1, par exemple:Generator_speed_std, Rotor_speed_std et Gearbox_bearing_1_temperature_std. De plus, nous avons quelques unes liées négativements corrélées, par exemple Nacelle_angle_std.
Egalement, nous avons des variables qui sont positivement et négativement avec les axes 2 et 3.
\(\bullet\) Contribution des variable à l’axe 1-2
Les variable avec une contribution supérieure à ce seuil représenté en pointillé rouge pourraient être considérées comme importantes pour contribuer à l’axe 1-2, par exemple Rotor_bearing_temperature, Gearbox_inlet_temperature et Rotor_speed.
Nous pouvons également utiliser le cercle de corrélation pour visualiser la contribution des variables à l’axe 1-2.
Le graphe ci-dessus illustre la contribution des variables à l’axe 1-2 selon la couleur, les variables représentées avec la couleur orange sont celle qui contribuent le plus à l’axe 1-2, par exemple Generator-bearing_2_temperature.
MAC_CODELa variable factorielle `MAC_CODE a quatre modalité (WT1,WT2,WT3,WT4), nous la recoder de sorte à avoir les modalités {1,2,3,4}.
#Régression et sélection de modèles
L’erreur de prédiction qui sera utilisée pour toute la suite sur le jeu de données test est définie comme suit :
Commençons par appliquer une régression linéaire à notre modèle Complet.
D’après le test de Fisher des hypothèses:
\(H_0:\beta_1=.........=\beta_{35}\) vs \(H_1: \exists j \in \{1,...,35\}/\beta_j \neq 0\)
Nous avons la \(p-value: < 2.2e-16\), nous rejetons donc l’hypothèse \(H_0\) à tout niveau usuel: parmi les 35 variables du modèle,il y’a celles qui sont significatives mais on ne sait pas lesquelles.
Nous pouvons construire un test de niveau \(\alpha=0.05\) de \(H_0\) vs \(H_1\) en combinant les 35 tests de student:
\(H_0(i):\beta_i=0\) vs \(H_1(i) \neq 0\).
La région de rejet combinée est de la forme : \(R=\bigcup_{i=1}^{35}R_i(\alpha_0)\).
Le test de student au niveau corrigé \(\alpha_0=\frac{0.05}{35}=0.001\).
D’après les sorties de la fonction summary(modèle) toutes les \(p_value < 0.001\), nous rejetons donc les hypothèse \(H_0(i) \forall \in \{1,...,35\}\) et toutes les variables sont significatives.
De plus nous avons \(R^2=R^2_{adj}=0.91\), alors \(90 \%\) de la variance expliquée de la variable TARGET est expliquée par ce modèle.
L’erreur de prédiction pour ce modèle est égale à 20785.07.
Dans le but d’avoir à la fois le modèle le plus petit possible et qui a une erreur de prédiction minimale, nous allons faire une sélection de variable exaustive en utilisant deux critères: \(C_p\) de Mallows et le \(BIC\).
Le critère `C_p$ de mallows est optimisé pour la ligne en haut du graphique, donc le modèle conservé est à 18 variables (plus la constante) suivant:
TARGET ~ Pitch_angle + Pitch_angle_std + Generator_speed + Generator_speed_std + Generator_bearing_1_temperature + Generator_bearing_2_temperature + Generator_bearing_2_temperature_std + Generator_stator_temperature + Generator_stator_temperature_std + Gearbox_bearing_1_temperature_std + Gearbox_bearing_2_temperature + Gearbox_bearing_2_temperature_std + Gearbox_inlet_temperature + Gearbox_oil_sump_temperature + Nacelle_temperature_std + Outdoor_temperature + Outdoor_temperature_std + Rotor_speed
Le critère \(C_p\) de Mallows est optimisé pour la ligne en haut du graphique, nous conservons donc un modèle à 18 variables (plus la constante).
Maintenant, après avoir retracé le plot de sélection de variable sur le nouveau modèle de 18 variables, elles seront toutes concervées.
L’erreur de prédiction correspondant au modèle seléctionné par le \(C_p\) de mallows est de 20190.95 et elle est inférieure à l’erreur de prédiction du modèle complet 20785.07. Le nouveau modèle explique mieux la variable TARGET.
Le critère BIC est optimisé pour la ligne en haut du graphique, nous conservons donc un modèle à 18 variables (plus la constante) et ce le même modèle retenue avec \(C_p\) de Mallows.
L’erreur de prédiction correspondant au modèle seléctionné par le BIC est de 20190.95 qui est sémilaire à celle obtenue avec le \(C_p\) de Mallows.
Donc pour la suite, nous allons travailler avec le modèle selectionné par les deux critères.
Dans cette partie, nous allons ajuster et prédire avec Ridge. Pour celà, il faut d’abord remplacer les prédicteurs qualitatifs par des “dummy” variables et créer un objet de type matrice contenant les prédicteurs à l’aide de la fonction model.matrix().
On applique la fonction cv.glmnet du package glmnet qui nous permet de lancer une cross validation sur un modèles RIDGE pour une grille de lambda. Cette fonction permet de sélectionner à partir des données une valeur de lambda qui minimise un critère de validation croisée.
En utilisant la fonction plot on obtient directement le graphique représentant des modèles et on extrait celui qui est associé au plus petit lambda.
Par conséquent, la valeur de lambda entraînant la plus petite erreur de validation croisée est 38.56808.
## 18 x 1 sparse Matrix of class "dgCMatrix"
## s0
## Pitch_angle 4.7089423
## Pitch_angle_std 4.6136524
## Generator_speed 0.2758431
## Generator_speed_std -0.3365382
## Generator_bearing_1_temperature -21.5934090
## Generator_bearing_2_temperature -2.8707309
## Generator_bearing_2_temperature_std -299.5077785
## Generator_stator_temperature 13.1159049
## Generator_stator_temperature_std 5.9459525
## Gearbox_bearing_1_temperature_std 122.6860940
## Gearbox_bearing_2_temperature 14.5425030
## Gearbox_bearing_2_temperature_std -168.0106894
## Gearbox_inlet_temperature -13.3329586
## Gearbox_oil_sump_temperature 1.7644989
## Nacelle_temperature_std -30.5889794
## Outdoor_temperature 2.7660319
## Outdoor_temperature_std 147.5024779
## Rotor_speed 29.1346919
Ci-dessous les coefficients \(\beta\) pour \(\lambda=38.56808\), aucun des coefficients n’est égale à zéro. Alors la régression de ridge ne permet pas une sélection de variables.
L’erreur de prédiction correspondant à la régression RIDGE est de 30947.94, elle est supérieure à l’erreur de prédiction du modèle complet 20785.07 et à celle obtenue par le modèle sélectionné par BIC et \(C_p\) de Mallows.
L’arbre présente un noeud qui se divise en deux branches (les branches 2 et 3), en identifant ceux qui ont des valeurs spécifques pour la variable Rotor_speed (noeuds de gauche nomée Rotor_speed ). Chaque noeud représente une question, la réponse non étant toujours dans la branche droite de l’arbre.
D’après la sortie de la fonction summary:
Le critère de split : Rotor_speed< 11.795 qui correspond au improve=0.7453816 le plus elevé fait qu’il soit le premier noeud de l’arbre.
Le nombre total d’instances pour le noeud : 319586.
La valeur prédite de la variable à prédire :172.91070.
La première table du summary (voir aussi cptable ci-dessous) montre comment l’erreur de prédiction (xerror) estimée par validation croisée évolue en fonction du paramètre \(C_p\), ici appelé complexity parameter (cp ou CP, mais rien à voir avec le Cp de Mallows). La colonne xstd donne une estimation de l’écart-type de l’erreur de prédiction. La colonne relative error évalue l’erreur d’ajustement.
Le summary décrit chaque noeud, on retrouve dans mean la valeur de la fonction de régression attachée à la feuille.
Plotcp() fournit une représentation graphique au résumé d’erreur validé de manière croisée. Les valeurs de cp sont comparées à la moyenne géométrique pour illustrer l’écart jusqu’à ce que la valeur minimale soit atteinte.
La meilleure erreur de prédiction est obtenue dans la ligne 8. Le cp optimal est 0.01.
| cp | nsplit | rel error | xerror | xstd | ||
|---|---|---|---|---|---|---|
| 0.0100000000 | 7.0000000000 | 0.0505808470 | 0.0506668353 | 0.0003283289 |
On affiche la valeur seuil du paramètre cp optimal 0.01.
Pour la valeur de cp seuil, on a le nombre de splits de l’arbre correspondant (nsplit), l’erreur croisée xerror et son écart-type xstd. On prend en général la première valeur de cp (i.e. la plus grande) qui est à moins de un écart-type du minimum de xerror. Une fois que la valeur de cp est choisie, on peut récupérer l’arbre correspondant par elagage de l’arbre avec le cp optimal 0.01.
Arbre simplifié : chaque noud représente une question, la réponse non étant toujours dans la branche droite de l’arbre. Chaque feuille est étiquetée par la décision associée et par l’effectif classe par classe des individus affectés à la feuille.
Sur le graphique ci-dessus qui fait une comparaison entre les valeurs prédites et les vraies valeurs, on voit que les points sont regroupés en lignes parallèles. En effet, à cause des capacités limitées de notre machine nous n’avons pas pu faire l’arbre maximale, ce qui fait qu’avec notre arbre de régression nous avons peu de feuilles qui ont de grands éffectifs (une valeur prédite peut correspondre à plusieurs valeurs reèlles). Donc l’arbre de regréssion que nous avons construit ne modélise pas bien notre variable à expliquer.
En conclusion, nous avons vu les étapes qui nous ont donné le meilleur résultat et mené à choisir un modèle. La description des données brutes nous a donné différentes idées. De plus, cela nous a permis de résoudre le problème causé par le dimensionnement de notre jeu de données. Cependant, la regréssion peut très certainement encore être améliorée en exploitant d’autres pistes,de même pour les résultats de l’arbre de régresion en construisant l’arbre maximale. En effet, on peut appliquer d’autres méthodes pour la sélection de variables, on peut également utiliser des algorithmes plus puissants (ex: RandomForest, XGBoost …etc) et jouer avec leurs paramètres afin d’avoir des modèles plus performants.